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Abstract 

We discuss how to characterize the behavior of a chaotic dynamical system 

depending on a parameter that varies periodically in time. In particular, 

we study the predictability time, the correlations and the mean responses, 

by defining a local-in-time version of these quantities. In systems where 

the time scale related to the time periodic variation of the parameter is much 

larger than the "internal" time scale, one has that the local quantities strongly 

depend on the phase of the cycle. In this case, the standard global quantities 

can give misleading information. 
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I. INTRODUCTION 



The forecast of the behavior of a system when its evolution law is known, is a problem 
with an obvious interest in many fields of scientific research. Roughly speaking, within this 
problem two main areas of investigation may be identified: 

A) The definition of the "predictability time". If one knows the initial state of a system, 

with a precision 5o = \5x{0)\, what is the maximum time Tp within which one is able 
to know the system future state with a given tolerance ^max ? 

B) The understanding of the relaxation properties. What is the relation between the mean 

response of a system to an external perturbation and the features of its unperturbed 
state Jll? Using the terminology of statistical mechanics, one wants to reduce "non 
equilibrium" properties, such as relaxation and responses, to "equilibrium" ones, such 
as correlation functions 0. 

A remarkable example of tjpe-A problem is the weather forecasting, where one has to 
estimate the maximum time for which the prediction is enough accurate. As an example of 
type--B problem, of geophysical interest, one can mention a volcanic eruption which induces a 
practically instantaneous change in the temperature. In this case it is relevant to understand 
how the difference between the atmospheric state after the eruption and the hypothetical 
unperturbed state without the eruption evolves in time. In practice one wants to understand 
how a system absorbs, on average, the perturbation (5/(t) of a certain quantity / - e.g. the 
temperature - just looking at the statistical features, as correlations, of / in the unperturbed 
regime. This is the so called fluctuation/relaxation problem 0. 

As far as problem-A is concerned, in the presence of deterministic chaos, a rather common 
situation, the distance between two initially close trajectories diverges exponentially: 

\dx{t)\ 5oexp{\t), (1) 

where A is the maximum Lyapunov exponent of the system From (|lD it follows: 
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T,~iln(^^^. (2) 

Since the dependence on 5max and 5o is very weak, Tp appears to be proportional to the 
inverse of the Lyapunov exponent. We stress however that @) is just a naive answer to the 
predictabihty problem, since it does not take into account the following relevant features of 
the chaotic systems: 

i) The Lyapunov exponent is a global quantity, i.e. it measures the average exponential 

rate of divergence of nearby trajectories. In general there are finite-time fluctuations of 
this rate, described by means of the so called effective Lyapunov exponent 7t(T). This 
quantity depends on both the time delay r and the time t at which the perturbation 
acted p. Therefore the predictability time Tp fluctuates, following the 7- variations 

i- 

ii) In systems with many degrees of freedom one has to understand how a perturbation 

grows and propagates trough the different degrees of freedom [0. In fact, one can be 
interested in the prediction on certain variables, e.g., those associated with large scales 
in weather forecasting, while the perturbations act on a different set of variables, e.g., 
those associated to small scales. 

iii) If one is interested into non-inflnitesimal perturbations, and the system possess many 

characteristic times, such as the eddy turn-over times in fully developed turbulence, 
then Tp is determined by the detailed mechanism of propagation of the perturbations 
through different degrees of freedom, due to nonlinear effects. In particular, Tp may 
have no relation with A @. 

In addition to these points, when the evolution of a system is ruled by a set of differential 
equations, 

dx/dt= f{x,t), (3) 

which depend periodically on time, 
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f{x,t + T) = fix,t) 



(4) 



one can have a kind of "seasonal" effect, for which the system shows an alternation, roughly 
periodic, of low and high predictability. This happens, for example, in the recently studied 
case of stochastic resonance in a chaotic deterministic system, where one observes a roughly 
periodic sequence of chaotic and regular evolution intervals 0. 

As far as problem-B is concerned, it is possible to show that in a chaotic system with an 
invariant measure P{x), there exists a relation between the mean response {6xj{T))p after 
a time r from a perturbation 6xi{0), and a suitable correlation function Namely, one 
has the following equation: 



where S{x) = — In P(a;). Equation ensures that the mean relaxation of the perturbed 
system is equal to some correlation of the unperturbed system. Of course since, in general, 
one does not know P{x), eq. provides only a qualitative information. 

In this paper we discuss how one has to reformulate the predictability problem, both of 
type A and of type B, for systems where the "seasonal" effects are important, i.e., when 
some relevant characteristic times change, in a substantial way, periodically in time. 

In section || we discuss how to approach the predictability problem in systems with time 
periodic evolution laws. Section |T| is devoted to the study of a toy model whose behavior, 
in spite of the simplicity of the model, catches the basic features of a system with seasonal 
effects. In section |^ we show the results of numerical experiments, that illustrate the 
relevance of the concepts introduced in section |T[ 

II. THE CHARACTERIZATION OF SYSTEMS WITH PERIODIC EFFECTS 

We consider systems in which one can identify two time-scales: the first, that we call Tg, 
can be thought as due to the coupling with an external time-dependent driver; the second, 
an "internal" time-scale Tj, characterizes the system in the limit of a constant external 



{SXj{T))p 




dS{x{0)) 
dxi 



p 



(5) 
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coupling, that we call the "stationary limit" . In the following, the external time dependence 
will be assumed periodic with period T, so that Te = T. 

If the system, in the stationary limit, is chaotic, one can take as the internal time-scale 
the Lyapunov time, i.e., the inverse of the maximum Lyapunov exponent, T/ ~ 1/A. 

In the case Te ^ Tj one can assume that the system is adiabatically driven through 
different dynamical regions, so that on observation times short with respect to the long 
external time, it evolves on a local (in time) attractor. If during its tour across the slowly 
changing phase space the system visits regions where the effective Lyapunov exponents are 
enough different, then the system shows up sensibly different predictability times, that may 
occur regularly when the driver is time-periodic. Consider, for instance, a slight modification 



of the Lorenz model [|TT|, which is the first geophysical dynamical system where deterministic 
chaos has been observed: 

dxjdt = 10 {y — x) 

dy/dt= —xz + r(t)x — y (6) 
dz/dt = xy — ^z 
where the control parameter has a periodic time variation: 

r(t) =ro- A cos(27rt/r). (7) 

Since this model describes the convection of a fluid heated from below between two layers 
whose temperature difference is proportional to the Rayleigh number r, the periodic varia- 
tions of r roughly mimic the seasonal changing on the solar heat inputs. When is close 
to the threshold r^r = 24.74, where in the standard Lorenz model a transition takes place 
from stable fixed points to a chaotic attractor, and the amplitude A of the periodic forcing 
is such that r(t) oscillates below and above r^r, a good approximation of the solution for 
very large T, may be given by 



x{t) = y{t) = ±\l^{rit) - 1) zit) = fit) - 1, (8) 

which is obtained from the fixed points of the standard Lorenz model by replacing r with 
r{t). The stability of this solution is a rather complicated issue, which depends on the values 



of To, A , and T. It is natural to expect that if Vo is larger than rcr the solution is unstable. 
In this case, for A large enough (at least Vo — A < r^r) one has a mechanism similar to 
stochastic resonance in bistable systems with random forcing 0. The value of T is crucial: 
for large T the systems behaves as follows. If 

T„ ^ nT/2 - T/4 (n = 1, 2, . . .) (9) 

are the times at which r{t) = Vcr, one can naively expect that for < t < Ti - when 
r(t) is smaller than - the system is stable and the trajectory is close to one of the two 
solutions (1), while for Ti < t < T2 - when r{t) > Vcr - both solutions (1) are unstable and 
the trajectory relaxes toward a sort of "adiabatic" chaotic attractor. The chaotic attractor 
smoothly changes at varying of r above the threshold Vcr, but if T is large enough this 
dependence in first approximation can be neglected. When r(t) becomes again smaller than 
Tcr, the "adiabatic" attractor disappears and, in general, the system is far from the stable 
solutions d^); but it relaxes toward them, being attractive. If the half-period is much larger 
than the relaxation time, in general, the system follows one of the two regular solutions 
for T2„+i < t < T2n+2- However, there is a small but non-zero probability that the system 
has no enough time to relax to and its evolution remains chaotic. Figure |^ shows the 
time evolution of the variable z, for To = 25.5 and A = 4, in the cases T = 300 (a) and 
T = 1600 (b). They provide an unambiguous numerical evidence that the jumps from the 
chaotic to the regular behavior (and the contrary) are well synchronized with r{t), with 
probability close to 1 when the forcing period T is very long (see fig. |l| b). On the other 
hand, for small value of T the system often does not perform the transition from the chaotic 
to the regular behavior, as seen in fig. |l] a. 

It is worth stressing that the system is chaotic. In both cases, in fact the first Lyapunov 
exponent is positive. 

It is rather clear from this example that the Lyapunov exponent is not able to characterize 
the above behavior, since it just refers to a very long time property of the system, i.e., a 
property involving times longer than T. A more useful and detailed information can be 
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obtained by computing a "local" average of the exponential rate of divergence for initially 
close trajectories. By this we mean an average which explicitly depends on the time to, 
modulus the external period T, to test the behavior of the system in the different states 
of the external driver. In this way one can make evident different behaviors, if any, of the 
system. 

We therefore define the mean effective Lyapunov exponent, for the time to G [0, T] and 
for a delay r, as 

1 N-l 

(7(r)),,= hm Y.^ito + kT,T), (10) 

where 

1 \5xit + r)\ 

is the local expansion rate, and 5x(t) evolves according to the linear equation 

From this definition it is clear that (7(T))t„ measures the growth of the distance between two 
trajectories that differ by when the external driver passes through a fixed value (or 

state). The maximum Lyapunov exponent of the system gives the global average of {'y{T))to'- 

A = lim (7(r)),„ = - / (7(r)),„dC (13) 



TJo 

If one is interested on predictability for times much smaller then T, {'y{T))t^ is a more 
appropriate quantity than A, since it distinguishes among different regimes. For example, in 
the system (H) discussed above, for the given values of the parameters, one has A > 0, but 
(7(r))t„ < when to G [{n — 1/4)T, {n + 1/4)T]. In the case of weather forecasting, different 
values of 7 for different to, correspond to different degree of predictability during the year. 

As far as the response properties are concerned, we expect that in a chaotic system 
the hypothesis of existence of "adiabatic" attractors implies that a fluctuation/relaxation 
relation holds also as a time-local property, provided one uses correlation and response 



functions computed according a local, not a global, average. So besides the usual correlation 
function between the variables Xi and Xj, 

C^f\r) = x,it)xj{t + T) -x-iX], (14) 

where (■) indicates the global average, 

A = lim - [' Mt')dt', (15) 



t^oo t Jo 

we introduce their correlation on a delay r after the time to G [0,T]: 

Cij{to,r) = {xi{to)xj{to + T))t, - {x^)t,{xj)t^, (16) 
where the local average is defined as in (p!OD 

1 N-l 

(A),„= hm -Y.Aito + kT). (17) 

In a similar way, one can consider two different kinds of mean response function of the 
variable Xj to a small perturbation 6xj: the global average response. 



and the local average response for the time to, 

The quantity ([T8|) gives the mean response, after a delay r, to a perturbation occurred at 
the time t, chosen at random, i.e., with uniform distribution in [0,T]. We shall see that 
R\j (r) can be rather different, even at a qualitative level, from Rij(to,T). 



III. A TOY MODEL 

The ideas discussed in Sec. || are illustrated here by means of a simple model. Consider 
the Langevin equation W2 
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j^q{t) = -a{t)q{t)+m, 



(20) 



where ^(t) is (5-correlated white noise, i.e., ^(t) is a gaussian variable with 



(21) 



and the coefficient a{t) is a periodic function of period T: a(t + T) = a(t). We require that 



dta{t) > 0, 



(22) 



to ensure a well defined asymptotic probabihty distribution for the stochastic process given 
by eq. (^). Moreover, we assume a slow variation of a(t), i.e. 

1 



mm a{t) > ^, 



(23) 



so that, by making an adiabatic approximation, a "local" probability distribution exists at 
any time. 

Without the noise term, the process described by eq. (^) is nonchaotic. Therefore the 



model (pOD cannot exhibit the full rich behaviour of chaotic systems, nevertheless it catches 
some of the relevant features. It is easy to see that the characteristic decay time of the local 
correlation 



Cito,T) = {q{to)q{to + r))t^ 



N-l 



lim - (ii^T + to + r) q{kT + Q 



(24) 



k=0 



depends on to. This can be easily computed by using the formal solution of 



q{t) = G{t) 



g(0)+ fdrG-\T)i{r) 
Jo 



(25) 



where 



G{t) = exp 



— / dr a{T) 





(26) 



A straightforward calculation leads to 



C{to, t) = C{to, 0) Gito, T)/G{to) (27) 
where the equal-time correlation is 



1 rt 

lim - Y: q{kTY + 2T rdrG-'ir) 



(28) 



In fig. H we show C(to, T)/C(to, 0) as a function of r for 



/27T 

a(t) = a + bcos[—t) , a + b>0 (29) 



T 



for two different values of to, namely to = and to = T/2, with T = 10, a = 1, 6 = —0.9 and 
r = 0.5. The different behaviour is evident. By defining a characteristic time as the time 
s it takes to have C(to,s) = 0.1, we get for this case Sq ~ 3.475 and st/2 ~ 1.275. When 
starting from to = T/2 the decay is almost a factor 3 faster than starting from = 0. The 
usual global average, 

C(^) (r) = lim - [' dt' q{t') q{t' + r) (30) 

t-»oo t Jo 

gives an average correlation function, so its characteristic decay time is not able to distinguish 
different regimes. Moreover, while C(to,T)/C{to,0) does not depend on the noise strength 
r, C(^)(r)/C(^)(0) does. In fig. we used T = 0.5. 

We consider now how the system responds at time to + t to a perturbation done at time 
to- This is described by the mean response function R{to,T), which can be computed as 
follows. One takes two trajectories differing at time to by a quantity e, i.e., 5q{kT + to)=e 
for any k, and evolving with the same realization of noise. Then the local response function 
is 

r,!, . y 1 5q{kT + to + T) 



N^^N ^0 5q{kT + to 



lim 1 E (31) 



k=0 

where bqikT + to + r) is the difference between the two trajectories at time to + t. Both 
times to and r run over a cycle, i.e., in the interval [0,T]. 
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By making use of (|25|) it is easy to see that 

RM^^. (32) 

By combining eq. (^T]) and (0) we have the fluctuations/relaxation relation p[ 

C{to,r) = C{to,0)Rito,T). (33) 

The scenario just described remains basically valid for non-linear Langevin equation. 
Consider for example the equation 

±q(t) = -a{t)q%t)+m (34) 



where a{t) is still given by 



It is natural to expect that, because of the adiabatic assumption, for eq. ( p4D a probability 
distribution 

Pt{q)(xexp-[Stiq)], (35) 

with St{q) = g''a(t)/4r, exists at any time. Therefore a natural ansatz is to assume that (^ 
becomes: 

In Fig. ^ we show i?(to, r) and ( g(to, t) [^St^{q)/^q(to)])^^ versus r, for different to, with 
a = 1, b = —0.5, T = 10 and F = 0.5. The two curves refer to to = and to = T/2. We see 
that eq. ( ^61) is well obeyed and that the characteristic decay times are different. 



IV. NUMERICAL SIMULATIONS 

In this section we report the results obtained for the Lorenz model (|^) with the Rayleigh 
parameter r varying periodically in time according to eq. (|^. We remind that at the value 
r = Tc = 166.07 the standard Lorenz model has a transition from a regular evolution (stable 
orbit) to a regime of intermittent chaos, and the maximum Lyapunov exponent depends on 
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r > Tc through a scahng law A ~ — Tc, for r — ^ 1 |T3[- Since we are interested in 



chaotic systems with a non neghgible variation of the degree of chaoticity, we chose Vo close 
to Vc- Let us recall that A changes very slowly for r near r^r = 24.74, so that a periodic 
variation of r in this region is not very interesting, at least if r{t) > for any t. We take 
To = 166.6, A = 0.5 and we consider three different periods: T = 1000, 100, 10. If T is very 
large with respect to the period (~ 0(1)) of the unstable periodic orbit, the variations of r 
can be taken as quasi-adiabatic. Therefore we expect that for large T and for t <^T, the 
local effective Lyapunov exponent of the system varies periodically in time with the same 
period of r,^e., 

~ A[r(t,)] ~ ^r{to) - r, oc ^ cos{2Tito/T). (37) 

From fig. | a one sees that this approximation holds true if r{to) is not too far from r^. 
The value of the maximum Lyapunov exponent should be close to the maximum Lyapunov 
exponent of the stationary Lorenz model computed at Tq, i.e., the average value of r on the 
whole period. 

When T is not long enough, so that one cannot speak anymore of adiabatic attractors, 
the behavior of (7)^0 has no relation with that of r(t), as shown in fig. ^ b. 



We discuss now the behaviour of the "local" correlation functions, defined by eq. (|16|), 
and of the "local" response functions, eq. (p^Of) . From fig. | one sees rather clearly how the 
"local" correlation functions depend rather strongly from the initial time t^. In particular, 
the dependence on to of the mean decay time is well evident. 

The correlation functions show a rapid oscillating behaviour on times ~ 0(1), with a 
typical period of the order of the mean circulation time near the "ghost" of the stable orbit 
existing at r = r^, while they decay for long delays, with characteristic times that sensibly 
depend on the degree of chaos relative to the initial instant. 

Let us stress that, because of the highly non gaussian nature of the system, there is 
not a simple direct proportionality between the response Rij and the standard correlation 
Cjj, both in the global and in the local version. This happens also in autonomous systems 
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10| and is not a pathological behavior. From very general arguments [|T0[ one can show 



that mean response and correlation function have the same qualitative behavior, e.g., the 
same decaying properties on large time delay. However, the agreement between Rij and Cij 
is very poor for moderate delay. This is so because in chaotic systems the error bars, in 
the numerical computation of the mean response function, increase exponentially with the 
delay and, when r is not very large, the mean response has to be compared with a suitable 
correlation function - see eq. (|^) - which depends on the unknown invariant probability 
distribution Pix). 



V. CONCLUSIONS 

We have discussed how to characterize the behaviour of a chaotic dynamical system 
when a "seasonal" effect is present, i.e. when there exist two well separated time scales: 
the internal one and that one of the periodic variation of a control parameter. A proper 
characterization has been obtained by the introduction of a restricted average for the rele- 
vant quantities: correlation functions, response functions and Lyapunov exponents. These 
selective averages keep into account, in an explicit way, the phase of the external period at 
which the system is observed. We stress that, in the presence of a "seasonal" effect, the 
usual global averaged quantities, can give only a rough information. 
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FIGURES 

FIG. 1. Model (|),(D with = 25.5 and A = 4. z as a function of t/T for (a) T = 300 and 
(b) T = 1600. 

FIG. 2. Model (M), (M) with a = l,b= -0.9, T = 10. and F = 0.5. C^'^\t) (diamonds) and 



C{to,T)/C(to,0) as functions of r for to = (full line) and to = T/2 = 5 (dashed line). 



FIG. 3. Model (H), (|2|) with a = I, b = -0.9, T = 10. and F = 0.5. i?(to,r) (full line) and 
{q{to,T) [dSta[q) / dq{to)]) (full squares) as functions of r, for to = and to = T/2 = 5. 

FIG. 4. Model (|),(|) with Vo = 166, 6 and A = 0, 5. {-i)t, as a function of to for (a) T = 1000 
and (b) T = 10. 

FIG. 5. Model (||),(0) with = 166, 6, ^ = 0, 5 and T = 1000. C33(to, r) as a function of r for 
(a) to = and (b) to = T/2 = 500. 
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